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Abstract 

We study chemically driven running droplets on a partially wetting solid substrate by means of coupled 
evolution equations for the thickness profile of the droplets and the density profile of an adsorbate layer. 
Two models are introduced corresponding to two qualitatively different types of experiments described 
in the literature. In both cases an adsorption or desorption reaction underneath the droplets induces a 
wettability gradient on the substrate and provides the driving force for droplet motion. The difference lies 
in the behavior of the substrate behind the droplet. In case I the substrate is irreversibly changed whereas 
in case II it recovers allowing for a periodic droplet movement (as long as the overall system stays far away 
from equilibrium). Both models allow for a non-saturated and a saturated regime of droplet movement 
depending on the ratio of the viscous and reactive time scales. In contrast to model I, model II allows for 
sitting drops at high reaction rate and zero diffusion along the substrate. The transition from running to 
sitting drops in model II occurs via a super- or subcritical drift-pitchfork bifurcation and may be strongly 
hysteretic implying a coexistence region of running and sitting drops. 
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I. INTRODUCTION 



The movement of droplets in external gradients fascinates scientists and layman alike at least 

since Newton's description yj] of Hauksbee's experiment with drops of orange oil that move be- 

J I 

tween two non-parallel glass plates towards the point of smallest plate distance [2J]. In another 
example a drop of liquid freely immersed in another liquid subject to a temperature gradient will 
move towards the higher temperature region due to Marangoni forces caused by surface tension 
gradients yfl. A drop sitting on a solid substrate also moves in a temperature Q4£| or wettability 
15, (1 7] gradient. Especially the Marangoni force is already used to manipulate droplets, for ex- 
ample in light induced drop movement [8Q. Similar concepts of a directed movement of small 
amounts of soft matter in a given gradient are also realized in models of cell motility [^]. 

However, even more intricate are situations where matter spontaneously starts a directed move- 
ment in initially homogeneous settings. Small pieces of camphor that move on a liquid surface 
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1 311 . More recently 



by emitting a surfactant have also been studied for centuries ifioi 111 , 
oil droplets containing volatile additives and interacting droplets of different volatile oils have 
been reported to move on solid surfaces due to the solutal Marangoni effect caused by evapo- 
ration/condensation [14]. Also intricate is the spontaneous movement of a juxtaposed pair of 



This ef- 



droplets with different wetting properties, so called bi-slugs, along a capillary tube 
feet was already mentioned by Marangoni llfjl . Another example are drops immersed in a second 
liquid. If the drops contain a soluble surfactant undergoing an isothermal chemical reaction at 



their surface the drops may start to move jln\ 



Apart from chemical reactions, drop movement 



can also be driven by surface phase transitions [|18L I9L The movement is possible because such 
active drops change their surrounding and produce a gradient that drives their motion. 

less wettable 

<Kx) 




substrate 



FIG. 1 : Sketch of a right moving droplet driven by a self -produced wettability gradient. 



Recent experiments also found chemically driven running droplets on solid substrates [2 
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22l,|23i|24j,|25L|260. In these cases small droplets of solution are put on partially wettable substrates. 
A droplet changes the substrate by adsorption or desorption of a solute rendering the substrate un- 
derneath the droplet less wettable than the bare substrate. The radial symmetry still assures an 
equilibrium position that is, however, unstable. As a result, fluctuations break the symmetry and 
the drop starts to move in a self-sustained manner. In the course of its movement it changes the 
substrate and leaves a less wettable trail behind (see Fig.Q]). We distinguish two types of experi- 
ments: In case I the substrate is changed in an irreversible way whereas in case II it recovers its 
initial state after the droplet has passed. The droplet does not rest again and follows a random 
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trajectory until trapping itself between self -produced less wettable patches (type I) [|21L 
ternatively, a running drop may perform a periodic movement (type II) 11241. 12511 . The droplet often 
moves with nearly constant velocity for a rather long time (tens of minutes). However, in all cases 
the motion ceases when the overall system reaches thermal equilibrium. Next, we shortly recall 
the specific results of different sets of experiments, yielding running droplets. 

In Ref. [|2l|l droplets of n-alkanes (n-octane and n-dodecane) are used that contain 1,1,2,2-tetra- 
hydroperfluorodecyltrichlorosilane (CF 3 (CF 2 )7(CH 2 )2SiCl3). The silane molecules form dense 
crafted monolayers on silicon or glass and render these surfaces hydrophobic. The droplet motion 
is stalled when no further hydrophilic surface is available. No limitation due to silane depletion in- 
side the droplet is observed. For millimeter sized drops (smaller than the capillary length) droplet 
velocities between 1 mm/s and lOcm/s are observed, depending on liquid viscosity, silane concen- 
tration and droplet size. The velocity increases with the silane concentration and the droplet size. 
The reaction rate is estimated to be 1 100 s _1 mol _1 . Typical reaction times a re g iven as 0.01-0.2 s. 

Another experiment is performed in a chemically different system [|22l I23Q using millimeter- 
sized droplets of a nonpolar solution of n-alkylamine (1 mM solution of C 6 NH 2 or Ci 8 NH 2 ) in 
decahydro-naphthalene (7 = 41mJ/m 2 , fj, = 0.001 Ns/m 2 ). Therein, silicon substrates with mi- 
croprinted high-energy surfaces are employed, that expose a dense packing of carboxylic acid 
functionalities (C0 2 H). The amines dissolved in the droplet adsorb at the substrate and produce 
a surface of lower energy exposing methyl groups. The effect of different adsorbates [22] and 
reaction kinetics ||23|l is investigated. The velocity decreases with increasing droplet size. 

The only example of a type II experiment [24] features droplets of oil (5 mM iodine solution of 
nitrobenzene saturated with potassium iodide) on glass substrates immersed in aqueous solution 
of a cationic surfactant (1 mM stearyl trimethyl ammonium chloride). In this system the stearyl 
trimethyl ammonium ion (STA) absorbs at the glass substrate outside the oil droplet and renders 
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the glass lyophilic. When the oil droplet moves on top of the coating the STA desorbs into the 
oil. In this way a wettability contrast between front and back of the moving droplet is created 
and sustained. In contrast to the type I case the substrate recovers its lyophilic state soon after the 
droplet has passed, because 'new' STA absorbs at the glass. The aqueous phase can be seen as an 
infinite reservoir of adsorbent. Movement only finally stops when the oil droplet is saturated with 
STA. 

Similar phenomena can be seen in metallurgic systems where droplets of liquid metals or alloys 
react with the metallic substrate, for instance, by alloying. The layer between the droplet and the 
substrate may be less wettable than the bare substrate resulting in the migration of reactive islands. 
This was studied for tin islands on copper surfaces that move and leave tracks of bronze behind 
1 2711 . Contrariwise, the layer can be more wettable than the bare substrate. This is typical for the 
related process of reactive spreading (also called reactive wetting) where a droplet of liquid on 



a (near' 
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y) non-wettable surface starts to spread after forming a more wettable layer underneath 



Mm 



33fl. Variants are possible in which during spreading the substrate becomes 



less wettable in the center of the drop j?>4 \. However, reactive spreading processes do normally 
not result in running droplets (but see the 'suddenly displaced' droplets in Ref. Q35U). 
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For type I experiments an implicit equation for the velocity v of the droplet was derived 11361 . 
from a simple theoretical argument. Based on a balance of friction force and driving capillary 
forces one obtains v = C tan#*(l — exp(—rL/v)), where r is the reaction rate, L the size of the 
droplet and C a constant. The dynamic contact angles at the advancing and receding ends of the 
droplet are then assumed to be identical (6*), i.e. the droplet profile is approximated by a spherical 
cap with 9* given by cos 9* = (cos 9 a e + cos 9 r e )/2. The static contact angles at the advancing edge 
9 a e and at the receding edge 9 r e > 9^ are different due to the chemical gradient. The expression for 
the velocity is found from a first order reaction on the substrate that yields chemical concentrations 
a a = and a r = 1 — exp(— rL/v) at the respective ends of the droplet. The expression for the 
velocity predicts a monotone increase of the droplet velocity with the droplet length L and the 
reaction rate r in line with experimental observations bill . The droplet velocity in the limi ting 
case of a saturated chemical reaction was also given in Q38[| . However, the experiments of Ref. Q23[| 
show the opposite trend; the velocity decreases with increasing drop sizes and effective reaction 
rate. In the framework of Ref. 113 611 the decrease is explained as a result of the flattening of the 
drops by gravity. Related works discuss running droplets in a random medium [3j3] and the forced 
wetting of a plate immersed into a reactive fluid [40] . 
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In this paper we propose and analyze dynamical models for self-propelled running droplets 
for both, type I and type II experiments. Our models consist of coupled evolution equations for 
the film thickness profile and the substrate coverage. Thereby the wettability of the substrate is 
modeled by a coverage dependent disjoining pressure. The different types of experiments are taken 
into account by incorporating different reaction terms. The models are capable of reproducing the 
different experimentally found regimes. A short account of a variant of the model for the type I 



experiments, i.e. for the case of a irreversibly changed substrate was recently presented 114 111 . There 
it was found that the dynamic contact angles at the advancing and receding ends of the droplet 
are not identical but rather resemble the corresponding static ones. Even for a strong driving 
force the deviation between dynamic and static angles is only about 10 percent. The bifurcation 
leading from sitting to running droplets for a finite reaction rate was identified as a drift-pitchfork 
bifurcation. Here we analyze the type I model in more detail and extend our analysis to the type II 
model. 

In section HI] the type I and type II dynamical models are presented. Results for the two models 
are discussed in section |In] and [IV] respectively. Stationary running droplets and sitting droplets 
are characterized in dependence of the control parameters reaction rate, droplet volume, diffusion 
constant of the adsorbate field and desorption rate of the adsorbate (for type II model only). The 
resulting families of solutions and their linear stability are used to derive phase diagrams that 
describe the existence regions for running and sitting droplets. To further elucidate the drift- 
pitchfork bifurcation that mediates the transition from sitting to running droplets also in the type II 
model we analyze the linear modes that destabilize sitting droplets close to the bifurcation. Before 
we conclude in section |Vj we use numerical simulations to illustrate and analyze periodic droplet 
motion that is possible in the type II model. The simulations are compared with our continuation 
results and the type II experiment of Ref. |24|]. 

II. MODEL 

A. Evolution equations 

The evolution of a thin liquid film on a horizontal smooth solid substrate is described by an 
equation for the thickness profile h(x, y, t) derived from the Navier-Stokes equations using long- 
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wave or lubrication approximation Il42h 



d t h = -V ■ l^V p\ . (1) 



u(x, z) = I — — zh ) Vp (2) 
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The parameters 7 and 77 are the surface tension and viscosity of the liquid, respectively. They 
define the viscous time scale r v = jL/r), where L is a typical length for the system. The change 
in time of the film thickness profile equals the gradient of a flow that results as the product of a 
mobility and a pressure gradient. The mobility h 3 /3r? corresponds to a parabolic velocity profile 

2 

within the film. The velocity component orthogonal to the substrate w can be obtained using the 
continuity equation 

d x u + d z w = . (3) 

The pressure 

p = 7A/1 + U(h, <p) (4) 
contains the curvature (or Laplace) pressure —'jAh and the disjoining pressure H(h,(j)). The 



latter comprises effective molecular interactions between the fi 



44, 



m surface and the substrate and 



45fl. As discussed in detail below, 



accounts for the wetting properties of the substrate ||43l . 
here a mathematically simple function H(h, </>) is used that is common in the literature. In many 
situations the qualitative outcome only depends on very general characteristics of the disjoining 



pressure Il46ll . The used form allows for solutions of Eq. (OQ) that represent static (i.e. sitting) 
droplets with a finite mesoscopic equilibrium contact angle. The disjoining pressure is chosen 
such that the droplets coexist with an ultrathin precursor film. 

The evolution of the density of the adsorbed layer on the substrate determines the wettability 
and is modeled by a reaction-diffusion equation for the dimensionless field < <p(x, y, t) < 1 

d t (f) = R{h,(f>) + d'A<f), (5) 

where the function R(h, a) describes adsorption or desorption on the substrate. The second term 
allows for a diffusion of the chemical species along the substrate. For simplicity we assume 
that the adsorbate has the same diffusion constant on the bare and the droplet covered substrate. 
However, assuming different constants would not change the presented results qualitatively. Note 
that for type I experiments </> directly corresponds to the substrate coverage. However, in a type II 
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experiment the droplet dissolves a more wettable coating. In this case the coverage corresponds 
to 1 — 0. Also here one can use Eq. ©, however, the signs of adsorption and desorption term are 
switched. With this convention, in both cases the wettability decreases with increasing 0. 

One can neglect the dynamics of the concentration field of the chemical in the bulk of the 
droplet by assuming a fast equilibration of the solute concentration within the moving droplet as 
compared to the reaction at the substrate. The fast equilibration is caused by diffusion and convec- 
tive motion within the droplet whereas the latter is driven by the lateral movement of the droplet 
along the substrate. It corresponds to the limit of a small Damkohler number Da= rcn/ (D/L) 
giving the ratio of reaction velocity at the substrate and diffusion velocity in the droplet 14711 • The 
parameter r is a typical reaction rate, c stands for a typical concentration of the chemical species 
in the droplet and D is the diffusion constant in the droplet. 



B. Reaction term 



Corresponding to the two different sets of experiments that the model shall describe we use two 
different reaction terms. For type I the initial coverage of the substrate is zero (i.e. = 0) and we 
only allow for adsorption underneath the droplet 

i*i(M) = nnt(h) (1-0). (6) 

The reaction saturates at = 1 because this is the maximal possible coverage. The function £(/t) 
represents a (smooth) step function that approaches one and zero inside and outside the droplet, 
respectively. For type II the 'rest' state of the substrate without any droplet is the fully covered 
substrate (again corresponding to = 0, see above section III A[) and we allow for desorption 
underneath the droplet and adsorption outside the droplet 

ifc(M) = r in £(h)(l-<f>) - tw[1 -£(/*)] 0. (7) 

The time scales of the reactions at the substrate inside and outside the droplet are defined by the 
effective rate constants r in and r out , respectively. Note that they have the dimension s -1 . The 
function £(h) may be the step function 

£i = Q(h-h B ), (8) 

or the smooth function 

6 = {tanh[(/ l -/i c )/A]+l}/2. (9) 
7 



The value of h c is chosen slightly larger than the thickness of the precursor film and the maximal 
drop height is always h max 3> h c . A small value of A h max ensures that the switch between a 
predominant adsorption reaction and a predominant desorption reaction occurs over a small film 
thickness range. Note, that £ 2 - > £1 f° r A — > 0. Changes in the details of the reaction term do not 
affect the results qualitatively. 



C. Disjoining pressure 



For the disjoining pressure we use throughout the present work 

2S«% , ( 0\ 5S s d 5 



^'^H 1 -v^ <10) 

where do = 0.158 nm is the Born repulsion length that defines a lower cut-off for the film thickness. 
The parameter Si and S s are the long- and short-range components of the total spreading coefficient 
S = Si + S s (for = 0). We use Si < and S s > corresponding to a destabilizing long-range 



van der Waals and a stabilizing short-range interaction [48]. For = the pressure allows for 
drops sitting on a stable precursor film as can be seen studying the corresponding densities of the 
excess surface energy /, related to the disjoining pressures by n = —d^f. 

The short-range interaction contains the influence of the coating and crucially influences the 



static contact angle [49]. To account for the varying wettability caused by the different substrate 
coverage we let the short-range part of the spreading coefficient S s depend linearly on the field 
(f)(x, y, t). The signs are chosen in a way that for both types of experiments g > assures that an 
increase in corresponds to a lower wettability, i.e. to a larger equilibrium contact angle 9 e given 
by cos# e = 5/7-1-1 |49|]. The constant g relates the coverage to the wettability and therefore 



defines the magnitude of the possible wettability gradient. Notejhat our Eq. (UOI) corresponds to 



the linear relation between cos 6 P and <t> assumed in Ref. 



D. Dimensionless equations 

We rewrite Eqs. ® to (flOl) by introducing scales 3777//K 2 , yfyTfH, and I for t, {x, y), and h, 
respectively. Then one obtains 2| S"; |c?§/Z 3 for the scaled spreading coefficient k. As length scale 
I we use the value of the film thickness where the local free energy f(h) has its minimum, i.e., 
where 11(h) = 0. This gives / = (55' s /2| ^ | ) 1 / 3 <io an d implies that the ratio of the strength of 
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the two antagonistic effective molecular interactions is only an implicit parameter of the system. 
The length / also corresponds to the thickness of the precursor film for a sitting drop on the bare 
substrate. 

Defining the dimensionless overall reaction rate r = 3r in 'jr]/lK 2 , diffusion constant d = 
3d'r]/Kl 2 , and ratio of reaction rates s = r in /r out , we obtain from Eqs. CQ) to (ITOT ) the dimen- 
sionless coupled evolution equations for the thickness profile h and the field 

d t h = -V {/i 3 V [Ah + U(h, </>)]} (11) 
d t cj) = rR(h,<f>) + dAcj) (12) 



with the disjoining pressure 



n(M) = -^+(i-£U <«) 



and the options for the reaction term 

iMM) = m (1-0) (i4) 

R 2 (h, 0) = £(/i) (1 - 0) - s [1 - £(/*)] 0- (15) 

To give an impression of the influence of the coverage we give in Fig. [2] droplet profiles for different 
values of the coverage 0. Thereby, is assumed to be independent of time and constant along the 
substrate. Increasing the ratio <f>/g clearly leads to increasing contact angles and decreasing droplet 
length. 

In the following we explore the full non-linear behavior of the dynamical models presented 
above. The main interest lies in an exploration of the existence regions of qualitatively different 
solutions depending on control parameters like reaction rate r, ratio of desorption and adsorp- 
tion rates s, and droplet volume V . It was already shown that a type I model allows for three- 



dimensional running droplets that move with constant velocity v and droplet shape 114 ill . Here, we 



restrict our attention to two-dimensional drops to be able to explore a large part of the parameter 



space. To determine sue 



ation techniques [5(1 51,1521 are employed to calculate two-dimensional running droplets moving 
with constant speed. This is achieved by switching to the comoving frame x — vt and imposing 
appropriate boundary conditions. In the comoving frame running droplets correspond to steady 
solutions. Integration in time of the full system is also used in some cases. Details on all numerical 
techniques used can be found in the Appendix. 



i droplets we therefore use Eqs. (II 11121) replacing V by d x . Then continu- 
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FIG. 2: Profiles of sitting droplets without reaction for different constant values of the wettability measure 
4>/g as given in the legend. The coverage is assumed to be independent of time and constant along the 
substrate. The system size is L = 10000 and all droplets have the volume V ~ 3 ■ 10 4 . 

III. RESULTS FOR MODEL I (ONLY ADSORPTION) 

We start with model type I, where adsorption takes place only underneath the droplet, i.e. 
we combine the evolution equations for the film thickness profile (fTT|) and the field of adsorbate 
coverage (fT2l) with the reaction term (fl4l) and the disjoining pressure (fl3l) . 

A. Thickness and coverage profiles 

First we focus on the characterization of the changing solution behavior in dependence of the 
overall reaction rate r. Without diffusion along the substrate (d = 0) one finds unstable sitting 
droplets and stable running droplets for all r. The droplets and the coverage profile move with 
constant velocity v and constant shape. We emphasize that this state corresponds to a subtle 
dynamical equilibrium given that coverage profile and droplet velocity depend on each other. 

Fig. [3] shows for two different reaction rates profiles of moving droplets [(a) and (b)] where the 
streamlines correspond to contour lines of the stream function 




(16) 
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and indicate the flow in the comoving coordinate system. Note that the velocity fields are obtained 
by (u, w) = (d z ip, —d x tjj). Also shown are the corresponding profiles for the coverage 4> [(c) and 
(d)]. 

The two sets of profiles belong to two qualitatively different regimes that are prominently visi- 
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es of the coverage (corresponding to the results obtained for a different disjoining 
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FIG. 3: Profiles of running droplets in the comoving frame. Shown are (a,b) the droplet profiles h and 
(c,d) the substrate coverage <fi for two different reaction rates r. Droplets on the left and right correspond 
to the non-saturated (r = 10~ 5 ) and saturated (r = 0.001) regime and move with velocities v sa 0.006 
and v rs 0.016, respectively. The streamlines plotted in (a) and (b) indicate the convective motion inside 
the droplets in the comoving frame and have the spacing Aip = 0.025 and 0.1, respectively. The direction 
of the movement is indicated by arrows. The remaining parameters are g = 2, d = 0.001, L = 10000, 
h = 3.8, h c = 2.0, Ah = 0.2 and the droplet volume is V « 3 • 10 4 . 



B. Dependence on reaction rate 

For low reaction rates [Figj3](a) and (c)] the coverage starts to increase at the advancing contact 
zone and continues to increase up to the receding contact zone where it is still well below the satu- 
ration value of <fi = 1. We call this the non-saturated regime. For high reaction rates [Figj3](b) and 
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(d)] the coverage starts to increase at the advancing contact zone as in the non-saturated regime. 
However, it increases much faster and reaches the saturation value 0=1 already underneath the 
droplet. At the receding contact zone it is always at the saturation value. We call this the saturated 
regime. The droplets in the two regimes behave qualitatively different when the reaction rate is 
changed as shown in Fig. SI 
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FIG. 4: Given is a characterization of stationary running droplets in dependence of the reaction rate r 
(logarithmic scale) for different strength' of the diffusion of the coverage 4> along the substrate d. Shown are 
(a) the velocity v for different diffusion constants d (see legend); and (b) the advancing (9 a , dashed lines) and 
the receding (9 r , solid lines) dynamic contact angles for d = 0.001 (heavy lines) and d = 0.1 (thin lines). 
The arrows with the labels r max , fdp and r sn in (a) indicate the locations for the maximum velocity, drift- 
pitchfork bifurcation and the saddle-node bifurcation, respectively, for d = 0.1. The remaining parameters 
are as in Fig.[3] 

For low reaction rates in the non-saturated regime the coverage does not reach the saturation 
value (f) = 1 at the rear of the droplet. This implies that the driving wettability gradient between 
front and rear of the droplet has not yet reached its maximally possible value. Therefore an increase 
in the reaction rate leads to a steeper increase in the spatial profile of the coverage underneath the 
droplet. In consequence a larger wettability gradient results implying a larger velocity. However, 
again it is subtle to determine the dynamic equilibrium. First, the higher velocity reduces the 
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contact time of the droplet with a given point of the substrate. Second, it also reduces the length of 
the droplet (see below). Both effects restrict the growth of the wettability gradient. The resulting 
increase of the velocity with increasing reaction rate can be seen in Fig.|4](a) for r < 0.0001. Note 
that it does only weakly depend on the diffusion along the substrate. 

Fig-Sl(b) shows the corresponding dynamical contact angles at the advancing and receding 
edges of the moving droplets. Below r « 0.0001 both angles increase, implying a decrease in the 
length of the droplet (that has constant volume). Because in the non-saturated regime the coverage 
in the contact zone at the advancing edge does effectively not depend on r, the increase in the 
advancing contact angle is caused by the larger velocity only. Seeing the full dynamic contact 
angle at a moving contact line as a superposition of a static (or equilibrium) part and a dynamic 
part, one can say that only the dynamic part of the advancing contact angle increases with r. 
The increase at the front is in line with the increase with velocity found for sliding drops on an 
incline (see for example Ref. J53]). However, at the receding edge the coverage in the contact 
zone depends strongly on r. Therefore the rather strong increase in the receding contact angle 
with increasing r is caused by changes in both, the static and the dynamic part. For a larger r the 
wettability at the back is smaller, i.e. the 'local equilibrium contact angle' increases. However, the 
dynamic receding contact angle normally decreases with increasing velocity [531]. This implies 
that the increase in the receding angle found here results from an increase of the static part only. 
It is even counteracted by a decrease of the dynamic part. Also in the contact angles, a relatively 
small d has no visible influence in the non-saturated regime. 

The maximum wettability gradient, i.e. driving force, and therefore the maximum droplet ve- 
locity is reached for the bare substrate ((f) = 0) at the advancing edge and maximum coverage 
(0 = 1) at the receding edge. The point of largest velocity [Fig.|4](a)] corresponds to the point of 
the largest difference between the receding and advancing contact angles [Fig.|4](b)]. We call the 
corresponding reaction rate r max . 

Next we discuss the behavior in the saturated regime (r > 0.0001 in Fig.[4|). For larger reaction 
rates the coverage at the rear of the droplet remains at its saturation value. However, the driving 
force and in consequence the velocity decrease slightly with increasing r (cf. also Ref. [|4l|]). Note, 
that for small d in Fig.|4](a) this is barely visible. The decrease from r rnax to r = 0.1 corresponds 
to about 2% of the maximal velocity. 

As detailed next this slight decay is caused by the dynamics in the advancing contact zone. 
In the saturated regime the time scale of the reaction is short compared to the one of the droplet 
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movement. Therefore the coverage increases very steeply underneath the droplet leading to an 
elevated coverage already in the advancing contact zone (see Fig.[3]b). The resulting decrease of 
the overall wettability gradient felt by the droplet causes the slight decrease of the velocity. In this 
regime the slight increase of the advancing contact angle results from an increase of the static part 
(because of the decreasing wettability) that overcomes a decrease of the dynamic part (decreasing 
velocity). The slight increase of the receding angle comes from the dynamic part only. 

We emphasize here that Fig.|4](b) clearly shows that the advancing contact angles are always 
smaller than the receding ones. As was already shown in Ref. |41] also here the differences be- 
tween the static and the dynamic contact angles at the front and the rear, respectively, are an order 
of magnitude smaller than the difference between the two static (or the two dynamic) contact 
angles. 



C. Phase diagrams 

1. Influence of Diffusion 



Let us finally discuss the influence of diffusion along the substrate on the dependencies shown 
in Fig.|H The result without diffusion (not shown) coincides nearly perfectly with the shown 
curves for d = 0.001. Only the decay in the saturated regime is slightly slower. We mention here 
that the decay also depends slightly on the used cutoff h c (cf. Eq. © and Ref. [41]). Increasing the 
diffusion has a rather small influence on the non-saturated regime but changes the saturated regime 
even qualitatively. The velocity decreases with the reaction rate because the adsorbate produced 
underneath the droplet close to the advancing edge diffuses onto the bare substrate ahead of the 
moving droplet. This effectively reduces the wettability there, implying a slower movement. In 
fact, running droplets cease to exist above a reaction rate r dp where the velocity drops to zero 
again, corresponding to a supercritical bifurcation for large diffusion (see curve for d = 1.0). 
The bifurcation structure can, however, be more involved. For a smaller d = 0.1 one finds that 
the branch of running droplets joins the branch of sitting drops in a subcritical bifurcation at 
r dp , i.e. decreasing v the running droplet branch first turns back in a saddle-node bifurcation at 
r sn and then joins the sitting droplet branch at r^. We indicate the definition of r sn and r<ip in 
Fig.|U(a) using the curve for d = 0.1 as example. However, the subcritical part of the branch is 
difficult to see. We refer the reader to the discussion of model II in section [IV] for more details. 
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FIG. 5: Existence of running and stationary droplets depending on the reaction rate r and the diffusion 
constant d. The solid line indicates a drift-pitchfork bifurcation where a running droplet solution branches 
off a stationary droplet solution. The dotted line indicates a saddle-node bifurcation where two running 
droplet solutions merge. Shown is also the transition between non-saturated and saturated regime (dashed 
line). The remaining parameters are as in Fig. [3] 

Between r dp and r sn both the large velocity running droplet and the sitting droplet are linearly 
stable and correspond to metastable states. The branch emerging from the subcritical bifurcation 
is linearly unstable and corresponds to threshold solutions that separate the two stable solutions. 
The prominent features (r max , ra p and r sn ) of the solution branches shown in Fig.@]can be used to 



determine phase diagrams by continuation techniques H52H . The phase diagrams show existence 
regions for the different types of solutions in the parameter space spanned, for instance, by the 
reaction rate and the diffusion constant (see Fig. [5]). To this end we followed the maximum of 
the v(r) dependence at r max , which marks the transition between the non-saturated and saturated 
regime for running droplets (dashed line in Fig. [5]). Continuation of the loci of the saddle-node 
bifurcation at r sn and the drift-pitchfork bifurcation at r dp gives the border of the existence region 
for stable running droplets (dotted line) and stable sitting droplets (solid line), respectively. The 
region between the two latter borders corresponds to a coexistence region where running and 
sitting droplets are metastable. Note, that unstable sitting drops do also exist everywhere in the 
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existence region of the running droplets. Such sitting droplets are steady solutions of the governing 
equations. However, because they are unstable, even infinitely small perturbations will grow with 
a characteristic rate. In consequence the droplets start to move and adopt the shape and speed of 
the stable running droplet solution. 



2. Influence of volume 

Practically, it is difficult to change the reaction rate over orders of magnitude. Experiments 



normally only cover a much smaller range fell 12311 . Nevertheless, the analysis of the v(r) de 



pendence (Fig.H]) for droplets of fixed volume gives a very good characterization of the general 
system behavior. An experimentally important parameter is the size or volume of the droplet. The 



velocity may increase 112 111 or decrease 112311 with increasing volume depending on the parameter 
regime 114 ill . We show in Fig. [6] the dependency of velocity on droplet volume (logarithmic scale) 
for a variety of reaction rates. Also there one can well distinguish a saturated and a non-saturated 
regime. In the non-saturated regime the droplet velocity increases with increasing size, whereas 
in the saturated regime it decreases. For small droplets one always finds a non-saturated regime 
whereas droplets of a very large size are always in a saturated regime (in Fig.|6]the latter is not yet 
reached for the curve with r = 10~ 5 ). The explanation for this behavior is closely related to the 
one given for the v(r) dependence. For small droplets the reaction does not reach saturation until 
the rear of the droplet has passed. This implies that an increase in droplet size gives more time for 
the reaction leading to a larger value of <p at the back and, in consequence, to a larger velocity. 

For very large droplets the reaction has enough time to reach saturation {<p = 1) at the rear, 
i.e. an increase in size does not increase the driving wettability gradient. It may even decrease 
slightly due to an increase of cf> in the contact zone at the advancing edge (see above). However, 
a larger droplet has a larger viscous dissipation leading to a decreasing velocity with increasing 
size. Inspecting Fig. [6] shows that for r > 5 • and V = 30 000 (i.e. the volume used in Fig.H]) 
one is well in the saturated regime. 

The volume V max where the maximal velocity is obtained can be followed in the space spanned 
by droplet volume and reaction rate. The obtained existence region for non-saturated and saturated 
running droplets are shown in Fig. [71 

Having studied the type I model describing experiments where the passing droplet irreversibly 
changes the substrate we next report on results for our type II model where the substrate recovers. 
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FIG. 6: The transition between non-saturated and saturated regime is most clearly seen in the dependence 
of the velocity v on drop size, i.e. droplet volume V (logarithmic scale). Curves are plotted for different 
reaction rates r (see legend). Remaining parameters are as in Fig.[3] 

In contrast to type I this allows for the description of a periodic droplet movement as experimen- 



tally observed in Ref. Il24ll 



IV. RESULTS FOR MODEL II (ADSORPTION AND DESORPTION) 

This section is devoted to the type II model that extends and generalizes the model of type I 
discussed above. Model II accounts for experimental situations where both, the droplet and its 
surrounding medium are able to change the substrate by adsorption or desorption. In this way a 
moving droplet makes the substrate less wettable, but after the droplet has passed the substrate 
may relax to its initial state. Such systems are modeled by extending the reaction kinetics for the 
0-field (fl4l) by an additional desorption term, i.e. we combine the evolution equations for the film 
thickness profile (fTTj) and the adsorbate field (IT2]) with the reaction term (fl3T) and the disjoining 
pressure (fT3l) . We will refer to this set of equations as type II model. The reaction term (fl"5l) is 
chosen such that adsorption and desorption of take predominantly place underneath and outside 
the droplet, respectively. 
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FIG. 7: Existence of running and stationary droplets depending on the reaction rate r and the droplet volume 
V (in units of 10 4 ) for two diffusion constants in the </>-field d = 0.001 (a) and d = 0.1 (b). The solid and 
dotted lines indicate a drift-pitchfork and a saddle-node bifurcation, respectively. The dashed line indicates 
the transition between non-saturated and saturated regime. Remaining parameters are as in Fig. [3] 

A. Thickness and coverage profiles 

As above for model I, we use continuation techniques to calculate running droplet solutions of 
model II as solutions stationary in a comoving frame. Examples of the resulting droplet profiles 
along with the profiles of the 0-field are shown for a very small diffusion d in Figs. [8] and [9j 
Thereby, Fig. [8] displays the results for several values of the reaction rate r for a fixed ratio s of the 
desorption to adsorption rates. The drop profiles change only slightly but the </>-profile undergoes 
prominent changes. As in model I there is practically no 0-field in front of the droplet (</> = 0). The 
field increases (i.e. the wettability decreases) underneath the droplet as in model I, but in contrast 
the 0-field decreases behind the droplet due to the desorption. The concentration of the 0-field 
at the rear of the droplet is increasing with r until it reaches saturation. Because r is the overall 
reaction rate also the desorption of cj) behind the droplet becomes faster, i.e. the tail of the 0-field 
becomes shorter. 

Fig. [9] displays droplet and profiles for a fixed reaction rate r but different values for the 
desorption to adsorption ratio s. The droplet shape and the form of the 0-field underneath the 
droplet are almost independent of s as one would expect. However the decay of the elevated 
4> value occurring behind the droplet is strongly influenced. For a small s, i.e. slow desorption 
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FIG. 8: Results of model II incorporating adsorption and desorption. Shown is the influence of the reaction 
rate r (see legend) for (a) the profiles of moving droplets and (b) the corresponding profiles in the (/>-field. 
The drop profiles at different r differ only slightly by visual inspection. The differences are more prominent 
in the ^-profiles. Parameters are s = 1, h c = 2, Ah = 0.2, g = 2, d = 0.001, V « 30000 (L = 10000, 
h = 3.8). 

as compared to adsorption, <\> decays slowly approaching qualitatively the behavior of model I. 
Increasing s reduces drastically the length of the <p tail behind the droplet. The length of the tail is 
very important when studying the periodic movement of a droplet on a finite stripe-like substrate 
(see below). 

B. Phase diagrams 

In the following the existence regions of running and sitting droplets are determined in their 
dependence on the control parameters reaction rate, desorption to adsorption ratio, diffusion con- 
stant and droplet volume. This is done along the lines explained in section Unl Continuation gives 
branches of stationary solutions depending on one control parameter. On these branches special 
points that separate qualitatively different behavior are identified and followed in the space of pa- 
rameters. The linear stability of the stationary solutions is also determined. Details on the used 
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FIG. 9: Shown is the influence of the ratio s of desorption and reaction (see legend) on (a) the profiles of 
moving droplets and (b) the corresponding profiles in the 0-field for model II. The drop profiles at different 
s differ only slightly by visual inspection. The differences are prominent in the tail of the 0-profile behind 
the droplet, r = 10 -4 and the remaining parameters are as in Fig. [8] 

techniques can be found in the Appendix. 

1. Influence of the desorption/adsorption ratio 

Studying first the influence of r and s we present in Fig.[TT](a) and (b) branches of stationary 
solutions in dependence of r for different values of s characterized by their velocity and the re- 
sulting existence regions of running and sitting droplets in the r — s parameter plane, respectively. 
Beside the moving droplets there exist sitting droplets (steady states) for all values of r and s. 
They are symmetrical with respect to the droplet maximum and may be stable or unstable. An 
example is shown in Fig.lTOl 

The curve for s — 1 in Fig.[TT](a) is used to better illustrate special values of the reaction rate 
introduced above in section Hill They mark the loci of the maximum of the v(r) dependence at 
r ma x, of the saddle-node bifurcation at r sn and of the drift-pitchfork bifurcation at ra p . As before, 
continuation in r of these loci when changing s gives the phase diagram Fig.[TT](b). In particu- 
lar one obtains the border between the non-saturated and saturated regime for running droplets 
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FIG. 10: Shown are profiles of (a) the film thickness and (b) the 0-field for running and sitting droplets (see 
legend) for model II. The drop profiles differ only slightly by visual inspection. The differences are more 
prominent in the 0-profiles. r = 10 -4 and the remaining parameters are as in Fig. [8] 

(r max ), the border of the existence region for stable running droplets (r sn ), and the border of the 
existence region for stable sitting droplets (r^). In the small region where both, running and sit- 
ting droplets are stable, one finds metastability. There noise, for instance, in the form of substrate 
inhomogeneities may lead to intermittent droplet movement. The sitting droplets existing in the 
gray shaded area of Fig.QTKb) are all unstable. 

Inspection of Fig.QT] reveals the rather strong influence of the desorption/adsorption ratio s 
especially on the transition from the non-saturated to the saturated regime and the slowing down 
in the saturated regime. The extension in r of the saturated regime is drastically reduced with 
increasing s. For decreasing s the maximum of the v(r) curve slowly transforms into a plateau and 
reaches for s = the form typical for model I. Comparing the strong decrease of the velocity with 
increasing r in the saturated regime [Fig.[TT|(a)] to the slight decrease found for s = (model I, 
Fig-H(a), curve for d = 0.001) poses the question why s has such a strong influence. 

The influence of s parallels the influence of the diffusion d in model I. There for large d the 
adsorbate diffuses in front of the advancing droplet rendering the substrate there less wettable. 
Thereby it decreases the overall wettability gradient between front and rear of the droplet. Here, 
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FIG. 11: (a) Velocity of running droplets depending on the reaction rate r. The branch of running droplets 
emerges from zero reaction rate and undergoes a subcritical drift pitchfork bifurcation at finite r > 0. 
The reaction rate-velocity curve is shown for three different values of the ratio between desorption and 
adsorption s. The labels r max , r^p and r sn indicate the locations for the maximum velocity, drift-pitchfork 
bifurcation and the saddle-node bifurcation, respectively, for s = 1.0. (b) Existence of running and sitting 
droplets depending on the reaction rate r and the ratio between desorption and adsorption s. Shown is also 
the boundary between non-saturated and saturated regime (dashed line). The remaining parameters are as 
in Fig.Hl 
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for large s the 0-field is removed behind the rear of the droplet rendering the substrate there more 
wettable. That also implies a reduction of the overall wettability gradient. 



2. The drift-pitchfork bifurcation 

To elucidate the mechanism of the transition between moving and sitting drops we focus for 
a moment on the drift-pitchfork bifurcation at r dp that separates the metastable and the running 



drop 



et region. 



This type of bifurcation is well known from reaction-diffusion (see, for instance, 



(1541 1551 . 1561 157L 15 811 and references therein) and hydrodynamic systems (see, for instance, 
60, |6jJ, [620 and references therein), where it mediates the transition between steady and travelling 
structures. In our system it breaks the reflectional symmetry of the sitting droplets leading to 
moving asymmetric droplets and acompanying travelling asymmetric adsorbate profiles. At the 
bifurcation a real eigenvalue switches sign, i.e. at r dp the velocity of the moving drops is zero. 
Beyond the bifurcation sitting droplets become unstable and start to move slowly. That makes it 
unlike a Hopf bifurcation (associated with the zero crossing of the real part of a pair of complex 
eigenvalues) where the travelling structure has a finite velocity at the bifurcation (for waves on 



flowing thin liquid films see the discussion in Ref. [63]). The bifurcation may be subcritical or 
supercritical (see Figs.|4l [TT1 and [T2l) and in consequence the bifurcating branch corresponds to 
unstable or stable moving drops. Approaching the bifurcation their respective velocity goes to 



zero as ^J\r — rd P \ providing an unambiguous signature of the drift-pitchfork bifurcation. 

The basic mechanism of the drift-pitchfork bifurcation is connected to the behaviour of the 
neutral (or Goldstone) mode related to the translational symmetry of the system. This mode with 
eigenvalue zero is obtained by analyzing the linear stability of the stationary solutions. In general, 
each continuous symmetry is related to such a neutral mode. In a sense, the neutral modes are the 
modes that are 'closest' to zero. This implies that a perturbation or modulation of these modes in 
an additional degree of freedom (if existing) gives modes that are probable candidates to cross zero 
and become instability modes. For instance, the transversal (fingering) instability of a liquid front 



results from a transversal modulation of the longitudinal translational neutral mode [16411 . Also 
one of the two coarsening modes of two liquid droplets is the combination of translational neutral 
modes of the individual droplets directed in opposite directions Q65Q . 

Here, the drift instability is associated with a mode representing a relative shift between the 
translational modes of the height and the <\> profiles. Right at the bifurcation this mode corresponds 
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exactly to the translational neutral mode. Although elsewhere one can still identify the translational 
mode when only looking at the sub-mode for the height profile or the one for the <p profile, looking 
at the complete mode one realizes that the relative weight of the two sub-modes is shifted in favor 
of one of them. This corresponds to the introduction of a relative shift between the two fields. The 
relative shift breaks the overall reflection symmetry and leads to the movement of the drops. 

3. Influence of diffusion 

Next we discuss the influence of r and d presenting in Fi g . [T"2l dependencies of droplet velocities 
on reaction rate for different diffusion constants and the resulting phase diagram in the r — d plane 
for a fixed desorption/adsorption ratio s. One notes first that the strength of diffusion has nearly no 
influence on the non-saturated regime and the transition value r max [the dashed line in Fig. [121(b) is 
practically vertical]. Note that Fig.[[2]focuses on relatively small diffusion. For large diffusion the 
droplets will also stop as shown for model I in Fig. [51 However, the saturated regime does depend 
on d quantitatively as well as qualitatively. For fixed r, with increasing d the stable running 
droplets become slower. In parallel its existence range in r shrinks slowly. As in model I this 
behavior is mainly caused by the diffusion of the 0-field in front of the advancing droplet. There 
it increases the coverage thereby reducing the overall wettability gradient, i.e. slowing down the 
droplets. The qualitative change concerns the character of the drift-pitchfork bifurcation. With 
increasing d it becomes less subcritical and at a critical d c it becomes supercritical, i.e. there 
exists no coexistence region for sitting and running drops any more. For small diffusion d < d c 
the existence range of stable sitting drops shrinks with increasing d whereas for larger diffusion 
d > d c it grows. 

4. Influence of volume 

Finally, in Fig.[T3lwe present the phase diagram for the dependence on reaction rate and droplet 
volume for fixed desorption/adsorption ratio s and diffusion constant d. The locations of the 
saddle-node, the drift-pitchfork bifurcation and the velocity maximum are all shifted slightly to- 
wards smaller r when increasing the droplet volume. The range in r of the non-saturated regime 
shrinks slightly, but the range of the saturated regime and the metastable region remain practically 
constant, they are only shifted towards smaller r. 
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FIG. 12: (a) Velocity of running droplets depending on the reaction rate r. The branch of running droplets 
emerges from zero reaction rate and undergoes a subcritical drift pitchfork bifurcation at finite r > O.The 
reaction rate-velocity curve is shown for three different values of the diffusion constant in the 0-field d. 
(b) Existence of running and stationary droplets depending on the reaction rate r and the diffusion constant 
of the (/>-field d. Shown is also the transition between non-saturated and saturated regime (dashed line). 
Remaining parameters are as in Fig. [8] 
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FIG. 13: Existence of running and stationary droplets depending on the reaction rate r and the droplet 
volume V (in units of 10 4 ). Shown is also the transition between non-saturated and saturated regime (dashed 
line). Remaining parameters are as in Fig.[8] 

Obviously the two contact regions and the substrate outside the droplet are not affected by a 
change in the droplet volume. However, an increase in droplet size increases the viscous forces 
and therefore stalls the droplet movement at smaller reaction rates. 



C. One-dimensional numerical simulations 



Finally, we employ numerical simulations of running droplets and show that the type II model 
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is able to describe different experimentally found modi of periodic droplet movement [24, 
To perform the simulations we use the routine 'd02cjc' provided by the NAG library Il66ll . It is 
based on a variable order, variable step size Adam's method. The simulations either utilize peri- 
odic boundary conditions (to model droplets on a ring-like track j24 \) or boundary conditions that 
mimick non-wettable borders of an otherwise wettable channel (to model droplets on finite stripe- 



like tracks) II24L I25H . Simulations are started from steady droplet solutions which develop in the 
absence of a chemical field, i.e. imposing = 0. Then the droplet movement is initiated by break- 
ing the symmetry of the 0-field by imposing a small gradient and starting the adsorption/desorption 
reaction. After a short initial transient the running droplets follow periodic trajectories that do not 
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depend on details of the initial symmetry breaking. For the periodic boundary condititions the ini- 
tial solution was one stationary droplet, whereas in the case of non-wettable boundaries the initial 
solution were two stationary droplets. 

In the case of periodic boundary conditions the droplets move with constant speed and shape 
after an initial phase. Fig.[14](a) shows space-time plots of the evolution of the film thickness 
for different reaction rates r. The droplet velocity increases with the reaction rate as expected 
from our continuation results. Fig.[TJ](b) shows a comparison of the droplet velocities obtained 
in the simulations and by continuation. The values for the simulations are estimated after the 
droplets have reached a constant speed and shape. One finds that both velocities match fairly well. 
However, for higher reaction rates the simulations slightly overestimate the velocities compared 
to the continuation results. This results from the lower and equidistant discretization used in the 
simulation in time. 



Experiments with droplets in a finite wettable channel found 'regular rhythmic motion' H24I1 or 
different types of 'shuttling motion' along with slowing and stopping behavior 12511 . The here per- 
formed simulations show a smooth transition between the different types of periodic movements 
depending on our control parameters. 

In general, the droplets in a wettable channel with non-wettable walls move periodically be- 
tween the two walls as shown for different reaction rates r in the space-time plots in Fig.[T5l(a). 
The initial solution of two stationary droplets quickly coarses upon starting the chemical reactions. 
The prevailing droplet oscillates between the channel walls with a frequency that depends on r. 

The droplet movement can be classified into two identical but antisymmetric halfcycles (i.e. 
with different signs in the velocity, and profiles that are related by reflection). Fig.[T5](b) shows 
the droplet velocity depending on time during one half cycle exemplary for two reaction rates r. 
We find that each half cycle typically contains three phases, distinguishable by their different 
velocities. After meeting the non-wettable boundary (phase I) the droplet velocity is very low or 
even zero in the case of very small reaction rates r. In this phase, the 0-field, which has been 
produced by thof e passing droplet has first to decrease, until the droplet can return on its own 
path. This phase of very small velocity is followed by a short phase of a very high velocity (phase 
II) until the droplet returns to a medium velocity (phase III), that is kept until it meets the opposite 
wall and the next half cycle begins. The subphases can also be very well distinguished in Fig.lT6l 
were we show hidden line plots of the film thickness profiles (top) and coverage profiles (bottom) 
for one period of droplet movement. One can clearly see that in phase I (velocity close to zero), 
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FIG. 14: Simulations of running droplets on a one-dimensional circular stripe (periodic boundary condi- 
tions). Shown are (a) space time plots of the film thickness for different reaction rates r (indicated on the 
right) and (b) a comparison of drop velocities obtained by numerical simulation (0) and continuation (solid 
line) for different values of the reaction rate r. Remaining parameters are g = 2, s = 1, d = 0.001, 
L = 100, h = 1.5, h c = 2 and Ah = 0.2. Dark (light) colors correspond to small (large) film thickness. 



after the droplet has encountered the non-wettable wall, the concentration in the 0-field is very 
high and drops steeply in phase II when the droplet starts moving again. Similar subphases of 
droplet motion have also been observed in Ref. [25] . In the context of that work the continuous 
transition from Fig.[l5](a) between r = 0.002 (top) and r = 0.0001 (bottom) can be seen as their 
transition between 'shuttling motion' and 'intermittent shuttling motion'. 



A simple bead-spring model based on a mechanical analogy is used in Ref. ||25|1 . It models the 
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FIG. 15: Simulations of running droplets on a one-dimensional finite stripe bounded by non-wettable bor- 
ders. Shown are (a) space time plots of the film thickness for different reaction rates r (indicated on the 
right) and (b) the velocity of the center of mass of the film thickness depending on time examplary for 
r = 0.001 (top) and r = 0.0001 (bottom). The thin dotted line in (b) at v = is meant to guide the eye. 
Remaining parameters are as in Fig. [14] Dark (light) colors correspond to small (large) film thickness. 

different experimentally observed regimes varying the wettability of the walls. Although it well 
captures the overall behaviour it can not resolve the more hydrodynamic aspects of the motion like 
the flow field inside the running droplets or the dynamical contact angles. 
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FIG. 16: Hidden-line space-time plots showing one halfcycle of the periodic movement of a running droplet 
on a one-dimensional stripe bounded by non-wettable borders. Shown are the film thickness (top) and the 
0-field (bottom), r = 0.0005 and remaining parameters are as in Fig.[l4l 

V. CONCLUSIONS 

In the present work we have developed and analysed two models for chemically-driven self- 
propelled running droplets on solid substrates. Such moving droplets were described for several 
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experimental systems using different liquids, substrates and reacting substances 



24, 



21, 



22, 



23, 



25Q. The movment of the droplets is driven by a self-produced wettability gradient that is 



perpetuated with the droplet itself by means of a desorption or adsorption reaction underneath the 
droplet. 



Two types of experimental systems were reported on: ( 



22, 



) adsorption underneath the droplet 



23Q; and (II) desorption underneath 



decreases the wettability of the substrate irreversibly J2I . 
the droplet removes a more wettable coating that is recovered behind the droplet through an ad- 
sorption from a surrounding medium 11241. 12511. We have described the droplet dynamics for both 
types of experimental systems using coupled evolution equations for the profiles of the film thick- 
ness and the adsorbate coverage. The equations can be derived from the Navier-Stokes equations 
using a long-wave or lubrication approximation [42] and assuming a small Damkohler number. 
The two types of experiments have been mapped onto two types of models. In the type I model a 
wettability-decreasing reaction takes place underneath the droplet. In the type II model this mech- 
anism is extended by additionally introducing a wettability-increasing reaction that takes place at 
the substrate outside the droplet. 

The wettability of the substrate enters both models through a disjoining pressure supplementing 
the Laplace pressure in the thin film equation. We have chosen here a disjoining pressure consist- 
ing of a long-range destabilizing part ~ h~ 3 and a short-range stabilizing part ~ /i -6 used, for 
instance, in Ref. [67] to study coarsening in dewetting. The long-range part corresponds to van der 
Waals interaction and is not influenced by the adsorbate. All the influence of the adsorbate goes 
into the short-range part that in the simplest case selected here depends linearly on the adsorbate 
coverage. 

Using continuation techniques and numerical simulations we have analyzed the solution be- 
haviour of both models. Thereby we have focused on stationary running and sitting droplets in 
two dimensions. Both models display a transition from a non- saturated to a saturated regime with 
increasing reaction rate. The transition is also obtained when increasing the droplet volume. In the 
non-saturated regime an increase in the reaction rate leads to a larger wettability gradient implying 
a larger droplet velocity. In the saturated regime an increase in the reaction rate does not increase 
the coverage at the rear of the droplet, i.e. it does neither lead to a larger wettability gradient nor 
to a larger droplet velocity. However, it has turned out that the driving force and in consequence 
the velocity are decreasing slightly with increasing reaction rate. This effect is due to a rise in 
the adsorbat concentration in the advancing contact zone. A similar behavior occurs when the 



31 



droplet volume increases. There, however, in the saturated regime the velocity clearly decreases 
with increasing volume because the constant driving force (wettability gradient) is counteracted 
by an increasing viscous friction. The latter dependencies found for the non-saturated and the 
saturated regime correspond very well to experimental results of Ref. 112 ill (Fig. 5) and Ref. Il23n 
[Fig. 7 (a)], respectively. To our knowledge there exist, however, no experimental results for a 
physico-chemical system that show the qualitative change between the two regimes in dependence 
of the droplet velocity on its volume. With the combination of materials used in Ref. [23] this 
should be possible because their Fig. 6 (a) shows the transition from the non-saturated to saturated 
regime for increasing solute concentration within the droplet. This corresponds directly to the 
dependence on the reaction rate shown here in Fig.|4j 

Allowing for diffusion of the adsorbate along the substrate does not affect the results in the 
non-saturated regime too much. However, it leads to a stronger decrease of the velocity in the 
saturated regime because adsorbate is transported to the substrate in front of the running droplet 
thereby decreasing the overall wettability gradient. This may even lead to a transition towards 
sitting droplets. This transition occurs either continuously through a supercritical drift-pitchfork 
bifurcation or discontinuously through a saddle-node bifurcation. In the latter case a metastable 
parameter region exist where running and sitting drops can coexist. 

We have determined the existence regions of running droplets depending on the reaction rate, 
the diffusion constant in the surface coating, the droplet volume and the desorption/adsorption 
ratio (only for model II). Both models show a strong dependence of the existence regions of sitting 
and running droplets on the diffusion constant and only a week dependence on the droplet volume. 
The strong dependence on the diffusion constant results mainly from the diffusion of the </>-field in 
front of the advancing droplet, thereby greatly affecting the driving force and the droplet velocity. 
The type II model also shows a strong dependence on the desorption/adsorption ratio, which is 
due to a decreased substrate coverage in close proximity behind the droplet, effectively decreasing 
the driving force. The primary effect of the desorption reaction is the possibility for the droplet to 
return on its own path, which we have illustrated by numerical simulations of droplets moving on a 
wettable stripe with non-wettable borders in an oscillatory manner, which has also been observed 
experimentally 11241 12511 . Finally, we have shown that different modes of periodic movement called 
'shuttling motion' and 'intermittent shuttling motion' in Ref. Il25n are covered by the presented 
model. 

The analysis of both models has shown that for the chemically-driven running droplets the ad- 
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vancing contact angle is always smaller than the receding one. As was already shown in Ref. [41] 
also here the differences between the static and the dynamic contact angles at the front and the rear 
are one order of magnitude smaller than the difference between the two dynamic contact angles. 
This challen ges the assumption of equal dynamic contact angles at the front and the rear that was 



used in Ref. [36, 



3711 to develop a simple description of self-propelled running droplets. A simple 
quantitative theory should instead be based on the assumption that the respective dynamic contact 
angles equal the different static contact angles at the front and rear. 

The proposed type I and type II model based on thin film or lubrication theory can very well 
reproduce the main features of droplet motion that have been observed experimentally. However, 
they fail to reproduce the reported damped oscillations in the droplet shape overlaying the contin- 
uous droplet movement Il25ll . These oscillations could be on the one hand the result of a weakly 
inhomogeneous surface, since the authors of Ref. [25] themselves suggested that they have no full 
control of the experimental surface properties. On the other hand the oscillations could be a sign 
of a Hopf bifurcation along the solution branch of running droplets. We are well aware that the 
models presented in this paper are minimal models that, however, reproduce qualitatively most as- 
pects of the behaviour of chemically self-propelled droplets. Refinements of the presented theory 
could include the viscous motion of the ambient medium as present in the type II experiments. 
This extension can still be based on the lubrication approximation, for instance, along the lines of 



the two-layer systems studied in Refs. [65, 68]. A second important extension could cover the case 
of a higher Damkohler number. Such a model has to include a description of the transport of the 
solute within the droplet. 

APPENDIX A: NUMERICAL TECHNIQUES 
1. Continuation 



Sitting and running droplets are steady solutions in the laboratory and comoving frame, respec- 
tively. They can be followed in parameter space using numerical continuation techniques 



for instance, employing the continuation software AUT097 Il52ll . The following section highlights 
some technical details of the employed techniques. 

The basic idea behind continuation is that unknown solutions of an algebraic system for a 
certain set of control parameters are obtained by iterative techniques from known solutions nearby 
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in parameter space. Differential equations of the form 

u'(x) = f(u(x),p) with f,u e R n (Al) 

subject to initial, boundary and integral constraints are discretized in space and then the resulting 
algebraic system is solved iterativly. Here the dash indicates the first derivative with respect to x 
and p denotes the set of control parameters. The presence of boundary conditions and/or integral 
conditions requires the presence of free parameters which are determined simultaneously and are 
part of the solution to the differential equation. The package AUT097 is limited to the continuation 
of ordinary differential equations (ODE's), thus it can only be used to compute droplet solutions 
in two dimensions. As an example we consider the continuation of stationary running droplets, 
steady in a comoving frame. After transforming Eqs. (fTT1) and (fT2l) into the comoving frame with 
velocity v and integrating the resulting time-independent thin-film equation we have the system of 
ODE's 

h[ = h 2 (A2) 
h' 2 = h 3 (A3) 

K = ^^-n.(/n,0i) (A4) 

<t>'i = 02 (A5) 
<f>' 2 = --(R^^ + vfc) , (A6) 

where /i is an integration constant. It has the physical meaning of a mean flow in the comoving 
frame. h\, h 2 , h 3 , 0i and 2 denote h, d x h, d xx h, <p and d x <p, respectively. The dashed quantities 
denote first derivatives with respect to x. The system is flux conservative, thus we need to specify 
the integral condition 

= - / h x dx-h (Al) 



L „ 

where L is the system length and h is the mean film thickness. Furthermore, for a complete 
description of the system we introduce periodic boundary conditions for the film thickness 

/n(0) = h x (L) (A8) 
h 2 (0) = h 2 (L) (A9) 
hs(0) = h(L) (A10) 
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and the mixed boundary conditions 







e(^ 1 (0))A 1 + {£(/*i(0)) + a [1 - ee»i(0))]} x 







(0 2 (O)-0i(O)Ai) 

e(^i(o))A 2 + + a [i - ecu w)]> x 



(All) 



(0 2 (L) - 0x(L)A 2 ) 



(A12) 



for the 0-field (see appendix I A 21) . Since translation of a running droplet solution in the x-direction 
also yields a valid solution one needs to introduce a pinning condition in the form of an additional 
boundary or integral condition, which will not be specified further in here. 

As mentioned earlier, the number of boundary conditions (NBCD) and integral conditions 
(NINT) imposes a constraint on the number of parameters (NPAR), that have to be varied dur- 
ing the continuation process. Specifically NPAR=NBCD+NrNT-NDIM+l, where NDIM is the 
dimensionality of the system of ODE's. This condition leaves us here with three free parameters, 
which are the principal continuation parameter (e.g. the reaction rate r), the mean flow // and the 
droplet velocity v. AUT097 uses the method of Orthogonal Collocation for discretizing solutions, 
where the solution is approximated by piecewise polynomials with 2-7 collocation points per 
mesh interval. The mesh is adaptive as to equidistribute the discretization error. Having specified 
the ODE system in standard form with boundary and integral conditions AUT097 then tries to 
find stationary solutions to the discretized system, by using a combination of Newton and Chord 
iterative methods. Once the solution has converged AUT097 proceeds along the solution branch 
by a small step in the parameter space defined by the free continuation parameters and restarts the 
iteration. 

The challenge usually is, to provide AUT097 with a nonuniform starting solution for the contin- 
uation. For our purpose it is sufficient to start the continuation close to the point of the primary 
bifurcation point, where the stable uniform solution of the ODE system (|A2I) - (1A6I) or a system 
similar to (IA2I) - (IA6I) with periodic boundary conditions undergoes a Hopf-bifurcation. In the 
vicinity of the bifurcation one can determine analytically small-amplitude sinusoidal stationary 
traveling waves. By selectively using reaction, boundary or integral conditions as primary con- 
tinuation parameter one finally computes fully nonlinear solutions for the film thickness and the 
coverage. 

AUT097 is not only able to follow solution branches but can also detect bifurcations, like saddle- 
node bifurcation or branching points, and can then follow these bifurcations in parameter space. 
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2. Boundary conditions 

The use of periodic boundary conditions does not rule out interactions between droplets in 
consecutive periods either via the film thickness or the 0-field. Interactions through the 0-field 
arise if the desorption of the coating is slow compared to the drop movement. One way to avoid 
this problem is to use very large periods, such that the 0-field has enough time to recover. This 
approach is successful in suppressing interactions but is very costly from a computational point 
of view. An alternative approach is to use other than periodic boundary conditions for the 0-field. 
The following paragraphs illustrate this approach. 

We generally consider a running droplet steady in a comoving frame, that has its maximum 
approximately in the center of the computational domain of length L. We assume that the film 
thickness profile obeys periodic boundary conditions hb = h(0) = h(L) and d z h(0) = d z h(L) <C 
1. The latter restriction ensures that the minimal film thickness is close to equilibrium and the 
period is large compared to the droplet. 

a. Model I 

In model I we assume that there is no chemical reaction taking place outside the droplet. There- 
fore, outside the computational domain for x < and x > L the following ordinary differential 
equation holds 

= d<f> xx + v(p x , (A13) 

which has the general solution 

0(x) = c + c'e - ^ (A14) 

with c and d being yet undetermined constants. For simplicity we assume that v > 0, i.e. the 
droplet is moving to the right. In front of the droplet it is assumed that — > as x — > oo. 
Furthermore we assume that the 0- field adopts a finite value behind the droplet such that — > 0^ 
as x — > — oo. First we consider the boundary at x = with the boundary value for the 0-field 
0(0) = 0o and the first derivative of the 0-field (f> x (0) = 4> x q. Since 0^ is a finite value we find 
d = 0, which leads to the Neumann boundary condition 

0,o = O. (A15) 
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Next we consider the boundary at x = L with the boundary value for the 0-field 4>(L) = 4>l and 
the first derivative of the 0-field <j> x (L) = 4> x l- Integrating ([A 131) we find 

= d<f) x + v<p + c" , (A16) 

where c" is an integration constant. Since <p — ► as x — > oo and also <fr x — ► as x — > oo we find 
c" = 0, which leads us to the mixed boundary condition 

= <f) L v + <j) xL d . (A17) 



b. Model II 



We assume for the /i-field outside the computational domain with x < or x > L h(x) = 
h(0) = h(L) — hh — const.. Then the concentration profile for the 0-field for x < or x > L 
can be computed analytically by solving a linear 2 nd order ordinary differential equation in x of 
the form 

= R 2 (h b , <p) + d(f) xx + v(f) x , (A18) 



which has the general solution 



where 



(x) = c + c'e XlX + c"e X2X , (A19) 



£(h b )+s[i-ah b )} 



(A20) 



r lv 2 r 



^ = -25 ± VlF + 5 { ^ ) + s|1 -« (ft ' )l} (A21) 

and d and c" are yet undetermined constants. For x — » ± oo <\> adopts a small but finite value 0oo. 

First we consider the boundary at x — 0. This leaves us with c" = 0. Using the boundary 
values for the 0-field 0(0) = O an d the first derivative of the 0-field X (O) = (p x0 at the boundary 
x = the following system of equations holds 

= fa-c-d (A22) 
= X . O - c% . (A23) 

Solving for d one finds the following mixed boundary condition 

= i{h h )X x + {i{h h ) + s [1 - ^h)}} (0,0 - 0oAi) . (A24) 
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The same procedure can be performed at the boundary x = L with <p(L) = <p L and 4> X (L) = 4> xL 
yielding the condition 

= £(h b )\ 2 + {£(h b ) +s[l- t(h b )}} {4> xL - L A 2 ) (A25) 

In the continuation algorithm the two periodic boundary conditions for the 0-field were substituted 
by the two mixed boundary conditions obtained above. 
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